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ABSTRACT 

Fourier series solutions to the deflection and magnification by a family of 
three-dimensional cusped two power-law ellipsoidal mass distributions are pre¬ 
sented. The cusped two power-law ellipsoidal mass distributions are characterized 
by inner and outer power-law radial indices and a break (or, transition) radius. 

The model family includes mass models mimicking Jaffe, Hernquist, and rj mod¬ 
els and dark matter halo profiles from numerical simulations. The Fourier series 
solutions for the cusped two power-law mass distributions are relatively simple, 
and allow a very fast calculation even for a chosen small fractional calculational 
error (e.g. ICC 5 ). These results will be particularly useful for studying lcnsed sys¬ 
tems which provide a number of accurate lensing constraints and for systematic 
analyses of large numbers of lenses. Subroutines employing these results for the 
two power-law model and the results by Chae, Khersonsky, & Turnshek for the 
generalized single power-law mass model are made publicly available. 

Subject headings: techniques: analytical — gravitational lensing — galaxies: 
structure 


1. Introduction 

Strong gravitational lensing effects provide us with direct tools to probe the potentials 
of lensing objects (i.e. galaxies and clusters/groups of galaxies) [see, e.g., Schneider, Ehlers, 
& Falco (1992) for an introduction]. Strong gravitational lensing effects can also be used 
to constrain cosmological models, in particular to directly measure the Hubble constant 
(Hq) in a standard cosmology via light travel time differences between lensed images of an 
extragalactic variable source (Refsdal 1964). Lensed image properties are related to (and 
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hence probe) the angular structures and radial behaviors of lensing mass distributions. 1 The 
angular structures and radial behaviors of lensing objects are usually studied using suitable 
parametric mass models. The choice of a mass model in lensing studies depends on two 
practical issues, apart from the question of whether the candidate model can be a good 
approximation to the “true” lensing mass distributions: (1) whether the parameters of the 
candidate model are constrainable with the available lensing data; (2) the availability of a 
method which can be used to calculate the deflection and magnification of the candidate 
model with the speed and accuracy required for the problems under consideration. 

This paper is concerned with the second issue for a family of three-dimensional triax- 
ial cusped two power-law mass distributions, which are described by p(R ) = poR~ Ui (l + 
R 2 )~(u 0 -vp /2 ^2 _ x 2 /a 2 + y 2 /b 2 + z 2 /(?) (Munoz, Kochanek, & Keeton 2001). This family 
of mass distributions include mass models resembling well-known galaxy models (e.g. Jaffe 
1983, Hernquist 1990, Dehnen 1993, Tremaine et al. 1994, Zhao 1996) and theoretical dark 
matter halo profiles from numerical simulations (e.g. Navarro, Frenk, & White 1997, Moore 
et al. 1998). These distributions allow more detailed studies of radial behaviors of galactic 
mass distributions via a break radius and an outer power-law radial index (i.e. v 0 in the 
above), compared to the generalized single power-law mass distributions (Chae, Khersonsky, 
& Turnshek 1998 [hereafter CKT]; Barkana 1998). It is also more realistic than the simpler 
two-dimensional double power-law models which were used by Evans & Wilkinson (1998) to 
study analytical properties of cusped lens models. To date, the only available method for 
calculating the deflection and magnification for this family of three-dimensional mass distri¬ 
butions has been based on numerically evaluating a set of one-dimensional integrals (Munoz 
et al. 2001; Keeton 2001). However, the numerical integration method is particularly slow 
for the cusped two power-law model family because the surface density for the model is not 
represented by a simple function (see §2). This paper presents mathematical solutions to 
the deflection and magnification by the cusped two power-law model family. The obtained 
mathematical solutions allow accurate calculations with errors of 10~ 6 or smaller for most of 
the parameter space and they improve the calculation speed dramatically for a large portion 
of the parameter space. 

This paper is organized as follows. In §2, we review the mathematical expressions for the 
projected surface mass densities of the three-dimensional cusped two power-law mass model 
and the two other models which are necessary to build up the model. In §3, the Fourier 
series formalism of CKT is reviewed and Fourier series coefficient functions of the three- 
dimensional cusped two power-law mass model are obtained for two parameter subspaces 


1 Gravitational lensing probes the total mass distribution of a lensing object, namely the sum of the 
luminous and dark mass distributions. 
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for which the radius is, respectively, smaller and larger than a break radius scaled with an 
ellipticity; mathematical details are outlined in the Appendices. A discussion of the results 
is given in §4. 


2. The Projected Surface Mass Density for the Cusped Two Power-law Model 


The three-dimensional cusped two power-law ellipsoidal mass model, referred to as ‘cusp’ 
model, is given by 


p cusp (X, Y, Z) = p 0 R~ Ui (l + R 2 )-^-^ 12 



( 1 ) 


where zy and zy, are inner and outer power-law radial indices respectively, and (A", Y, Z) 
and (a, 6, c) are body coordinates and “break sizes” along the axes, respectively. 


The surface mass density for the ‘cusp’ mass model projected on the lens plane is, up 
to a multiplicative constant controlling the deflection scale of the lens, 2 given by 
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where 

C = —[1 + e cos 2(0 — 0 O )] 1/2 , (3) 

Tb 

where parameter e(> 0) is related to the ellipticity (e) via e = 1 — [(1 — e)/(l + e)] 1 / 2 . In 
equations (2a) and (2b), B(a, b) is the beta function and F(a,b]C]z ) is the hypergeometric 
function. The relations between the parameters of the surface mass density [eqs. (2a) & 
(2b)] and the parameters of the three-dimensional mass density [eq. (1)] and three Eulerian 
angles can be obtained, e.g., following CKT. The Erst expression given by equation (2a) 
is suitable for calculating the deflection and magnification for r > rb/{ 1 + e) 1 / 2 , while the 
second expression given by equation (2b) is suitable for r < /(1 + e) 1 / 2 . 

Each term in (the expansion of) the right-hand side of equation (2a) is the surface mass 
density of the familiar softened general single power-law mass model (CKT; Barkana 1998) 
given by 

K - te (C;ri= (1 + ) y i , (4) 

for which Fourier series solutions for the deflection and magnification are given in CKT. As 
indicated by a superscript, the surface density of equation (4) is referred to as ‘sple’ (i.e. 
softened power-law ellipsoid) model. Equation (2b) consists of terms of the following form, 
which will be referred to as ‘Nuker’ model following the Hubble Space Telescope group who 
observed the centers of nearby early-type galaxies (e.g. Lauer et al. 1995), although the inner 
power-law index [2A in eq. (5)] can take positive values here, 

a2A 

^ Nuker (C;A,h)= (1+ c c2) , +1 , (5) 

where A > — 1. Note that multiplicative constants were suppressed in both equations (4) 
and (5) for mathematical convenience, as in equations (2a) and (2b). For positive integer 
values of A, equation (5) can be written as 

K N “ k “'(C; A,//) = £(-!)' ( )V p1 '(C;a<-A + /). (6) 

1=0 ' ' 

Summary of the mass models given above can be found in Table 1. 


3. Fourier Series Coefficients for ‘CUSP’ Model 

As was noticed in CKT for the ‘sple’ model [eq. (4)], the surface potential, deflection and 
magnification for a surface mass density with elliptical symmetry [i.e., which is a function 
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of C given by equation (3)] can be written as Fourier series with coefficients which can be 
written as linear combinations of the following “/-functions” 3 of r whose general definitions 
are given in Appendix A: /(°)(r), /(°')( r ) ; 4 42 (r), 42(u)> and /^(r) (m = 1,2,3, ...). 5 
Using these /-functions, we may write the surface potential, ip = ip(r, <f>), and its first and 
second derivatives in polar coordinates as follows: 


dip 
dr 


dip 
rdcj) 


d 

/ dip 

dr 

\ dr 

d 

( dip 

rdcf) 

\ dr 

I! 

' dip 

dr \ 

j'd(f) 

^ f 

' dip 

rdcf) \ 

K rd<p 



52 ^[42(4 + 42(0] cos[2m(0 - 0 O )], 

m=l 


J (0) ( r ) + ^[/W(r) - 42CO] cos[2m(0 - 0 O )], 

m=l 


52t42('0 + 42(0] sin[2m(0 - fa)], 

m= 1 


—/(°)(r) + / (Q,) (r) 


r 


(0 

( 8 ) 

(9) 


-- 52t( 2m + o42o) + (2m - o42(0 - 442(0] cos[2m(</> - </> 0 )], 

ra=l 

( 10 ) 

0| 52 m [42('0 - 42(4] sin[2m(0 - <&,)], (11) 

ra=l 
^ oo 

-- + 1)42(0 - (2m - 1)42(0] sin[2m(0 - 0 O )], (12) 

ra=l 
2 °° 

- 52 m [42('0 + 4m00] cos[2m(0 - 0o)]- (13) 

m= 1 


As shown above, the m-th order Fourier coefficients of the potential, its first and second 
derivatives are entirely determined by the m-th order /-functions. For any surface mass 
density with elliptical symmetry and a reasonable ellipticity, the magnitude of an /-function 
decreases rapidly as order m increases. The I^\r) and I^ 2 \r) functions are determined 
solely by the mass distributions within and outside the radius r, respectively. The /®(r) 
functions arise in the derivatives of the /^(r) and I^ 2 \r) functions. 

Below we obtain mathematical expressions for the functions /^(r), A°4(r), 42(r), 


3 These /-functions are constants for singular isothermal mass distributions. 

4 This symbol is used throughout in place of used in CKT. 

5 For to = 0, /(°l(r) = 42 ( r ) an( l A 0 \r) = 242 ( r )- 
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/q\ /ot 

1 2 m \ r ), and ^ 2 m \ r ) °f the ‘cusp’ model using the two expressions of the surface density 
[equations (2a) & (2b)]. As was mentioned in §2, the surface density expressions of the 
‘cusp’ model consist of those of the ‘sple’ and ‘Nnker’ models. Thus, in order to calculate 
the /-functions of ‘cusp’ model, we need the /-functions of the ‘sple’ and ‘Nuker’ models. 
Mathematical details outlining the derivations of the necessary /-functions of the ‘sple’ and 
‘Nuker’ models are given in the Appendices A and B. 


3.1. The Case r < r b /( 1 + e) 1 / 2 
From equations (2b) and (5), we have 
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and /^)) ktl ^'(r; A,/.) are the /-functions of the ‘Nuker’ model [eq. (5)]. In equations (15) 
and (16) (a ) b is the pochhammer symbol, i.e., (a) b = T(a + b)/T(a). The /-functions of the 
‘Nuker’ model are given below. 

a. A = non-integer 

As outlined in Appendix A, the /-functions of the ‘Nuker’ model for non-integer values 
of A are given by the following simple expressions, 
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where 
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b. A = zero or positive integer 

For zero or positive integer values of A, the /-functions can be calculated using equation 
(6) and the /-functions of the splc model, i.e. 
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where /^^(r; /i) are the /-functions for the ‘sple’ model given by equation (4). 


3.2. The Case r > r&/(1 + e) 1 / 2 

For this case, the /-functions are calculated using equation (2a), which may be re-written 
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where / 2 s ^J e ^(r; /i) are the /-functions of the ‘sple’ model [eq. (4)] for r > r&/( 1 + e) 1//2 which 
are given in CKT. The series sum of equation (26) can be efficiently calculated by substituting 
the series expressions for fi) given in CKT into the equation, and turning double 

infinite sums into single infinite sums using the following rule, 

OO OO OO j 
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where j — i + k. 

The following equations (28), (29), (30b), (31), and (32) are obtained from equations 
(26b), (43), (27b), (28), and (44) of CKT, respectively using the rule given by equation (27), 
while equation (30a) is obtained from a new expression for /^^(c n) which is given in 
Appendix B [eq. (B2)]: 
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In equations (28), (29), (30b), (31), and (32), the functions h(r), £i(r), and e 2 (r) are given 
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In equations (28) & (30b), the constant infinite sum can be evaluated using F(a, 6; c; 1) = 
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where m is zero or positive integer. Equation (30a) converges in a non-trivial fashion and 
requires a careful treatment. The convergence property of equation (30a) and a strategy for 
its efficient numerical evaluation are discussed in Appendix C. 


4. Discussion 

We have obtained Fourier series solutions to the deflection and magnification by the 
three-dimensional cusped two power-law ellipsoidal mass (‘cusp’) model [eq. (1) or eqs. (2a,b); 
see also Munoz et al. (2001)] by calculating the /-functions of the model whose general 
definitions are given in Appendix A. The /-functions of the ‘cusp’ model given in §3 are well- 
defined rapidly converging series and can be evaluated efficiently using recursion relations 
for gamma functions and those for hypergeometric functions for the successive terms in the 
series (Abramowitz & Stegun 1964; Gradshteyn & Ryzhik 1994). Using these /-functions, the 
deflection and magnification can be calculated quickly and accurately for arbitrary parameter 
values. In particular, for any value of the radius sufficiently smaller or larger than the break 
radius the calculation is very fast, even for highly accurate calculation, e.g., with fractional 
errors of 10 -5 or smaller. For example, compared with that for the ‘sple’ model for a given 
value of the normalized radius r/'/y, the calculational speed for the ‘cusp’ model is slower 
by factors of ~ 3 and ~ 8 respectively for r/r b = 10 and 0.1. 6 It typically takes ~ 10~ 5 
to ~ 10“ 3 seconds on an AMD Athlon 1.1 GHz CPU and an UltraSparc II 400 MHz CPU 
to calculate the deflection and magnification with fractional calculational errors of 10 -5 or 
smaller for any ellipticity e < 0.7. 

Based on the present results, it is only around the break radius (specifically for 0.7r fe /vT+~e < 
r < 1.3r'ft/\/l — e) that the calculation is limited, mainly due to difficulty in efficiently evalu¬ 
ating an / ( ^-function given by equation (30a). There are two limitations: First, for elliptici- 
ties larger than e ~ 0.7, smallest possible fractional calculational errors are rs-/ 10 5 . Second, 
the calculation slows down significantly compared with calculations outside the above range 
of radius; on the same machines and with the same limits on the fractional calculational error 
and ellipticity given above execution time is typically ~ 0.5 x 10~ 3 to ~ 10~ 2 seconds. These 
limitations, although not very significant, can be reduced if a new /(^-function simpler than 
equation (30a) is obtained in the future. 

To summarize, the present Fourier series solutions to the deflection and magnification 


6 The calculational speed in the Fourier series approach (for any model) depends on the mass ellipticity 
for any value of r/rb . Execution time increases approximately linearly with ellipticity for a fixed series 
truncation criterion. 
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by the ‘cusp’ model improve the calculational speed by factors of tens for a large portion 
of the parameter space of the model, compared to numerical integrations. For radii close 
to the break radius, the Fourier method still appears to be a few times faster than the 
numerical integration method for relatively small ellipticities (e.g. e < 0.5). The present 
Fourier series solutions for the ’cusp’ model can also be used to calculate very accurate 
the deflection and magnification with fractional calculational errors of 10 -6 or smaller, for 
the entire parameter space but the small portion of the parameter space specified above. 
Therefore, for practical applications of the ‘cusp’ model in which one is mainly interested in 
calculating the deflection and magnification at radii not too close to the break radius, the 
present Fourier series solutions greatly reduce his/her numerical efforts. 

Numerical implementation of the /-functions of the ‘cusp’ model given in this paper 
is relatively straightforward except for two of them which are discussed in Appendix C. 
However, for the convenience of the lensing community well-tested codes both for the ‘cusp’ 
[eqs. (2a,b)] and ‘sple’ [eq. (4)] models are made publicly available on the following web page 
http://www.jb.man.ac.uk/~chae. The codes for the ‘cusp’ and ‘sple’ models are given 
the names of ‘cuspFS’ and ‘spleFS’ respectively, where ‘FS’ stands for Fourier series. 

The remarkable calculational speed for a chosen small calculational error and the easy 
control over the calculational error (or, speed), which the Fourier series solutions both to 
the ‘sple’ and ‘cusp’ models demonstrate, suggest that the Fourier series method is the 
most efficient method for calculating the deflection and magnification for a general class of 
mass distributions. In particular, for any surface mass density which is a function of </ [eq. 
(3)], one only needs to calculate the /-functions of the model [eqs. (Al) through (A5)] to 
calculate both the deflection and magnification. As shown in CKT and in this paper, the 
/-functions both of the ‘sple’ and ‘cusp’ models are well-defined, relatively simple converging 
series. It would be expected that for other classes of mass distributions which have not yet 
been explored, similar relatively simple expressions for the /-functions exist. For example, a 
generalized softened power-law ellipsoid with an exponentially declining envelope given by 

m = (1 + % )y/2 exp (~R 2 /R: ! t ), (37) 

where R 2 = ( X/a ) 2 + ( Y/b ) 2 + ( Z/c ) 2 , has the projected surface mass density given by the 
following equation: 

<0 = «»(! + C 2 )' 1 - 172 exp (—(*/%)* (i,| ~ , (38) 

where T(a, 6; z) is a confluent hypergeometric function and ( is given by equation (3). Equa¬ 
tion (38) can be expanded in terms of ‘sple’ models, and thus the /-functions of the model 
can be calculated in a straightforward manner. 



- 12 - 


Before we conclude this paper, below we discuss astrophysical motivations and aspects of 
using the ‘cusp’ double power-law model family as lens models, which have not been empha¬ 
sized in previous work on gravitational lensing. For the ‘cusp’ model given by equation (1), 
the break radius and the outer power-law index have been used as defining outer envelopes 
of systems of stellar objects (e.g. Jaffe 1983, Hernquist 1990, Dehnen 1993, Tremaine et al. 
1994, Zhao 1996, Munoz et al. 2001), v a = 4 being a preferred choice for the outer power-law 
index. Thus, the break radius would be a scale size of the system and the inner power-law 
index would represent an overall radial profile within the scale size. However, we also would 
like to use the break radius as a means to allow the radial profile of the lens to vary from 
a shallower inner profile, which is on average shallower than isothermal, to a steeper outer 
profile which is on average similar to isothermal and not steeper than v Q = 3 in general. This 
would be supported both by observational evidence and theoretical arguments. First, rota¬ 
tion curves of spiral galaxies, dwarf galaxies, and low surface brightness galaxies generally 
rise in the inner regions of the galaxies and are nearly flat beyond a turnover radius (e.g. 
Kravtsov et al. 1998; Rubin, Waterman, & Kenney 1999; Sofue et al. 1999; van den Bosch et 
al. 2000; S waters, Madore, & Trewhella 2000; de Blok et al. 2001). Rotation curves of small 
samples of cllipticals/SO galaxies and blue compact galaxies show similar radial behaviors 
(Bertola & Capaccioli 1975, 1977, 1978, Rubin, Peterson, & Ford 1980; Ostlin et al. 1999). 
These rotation curves imply shallower-than-isothermal inner profiles and steeper outer pro¬ 
files which are similar to isothermal but generally not steeper than v Q = 3. Second, Hubble 
Space Telescope luminosity profiles of the centers of early-type galaxies show shallower-than- 
isothermal profiles (e.g. Rest et al. 2001; Faber et al. 1997; Gebhardt et al. 1996) virtually 
for all the galaxies imaged. These observations imply that, at least at the very centers of the 
galaxies where dark matter contribution to the total mass density is minimal, mass profiles 
are shallower than isothermal. Since we do not expect such shallow profiles to extend to 
large radii, varying radial profiles with relatively small break radii appear to be inevitable. 
Third, numerical simulations of halo formation in cold dark matter models predict mass dis¬ 
tributions with varying radial profiles (e.g. Navarro, Frenk, White 1997; Moore et al. 1998; 
Jing & Suto 2000; Ghigna et al. 2000) with shallower-than-isothermal inner profiles and an 
outer profile of v Q = 3. Moreover, mass distributions in the inner regions of the galaxies are 
expected to be modified due to the effects of the baryons (e.g. Blumenthal et al. 1986; Mo, 
Mao, & White 1998; Cole et al. 2000; Gonzalez et al. 2000; Kochanek & White 2001). It is 
thus of considerable interest to use gravitational lenses to determine the inner profiles, break 
radii, and outer profiles (and eventually truncation radii, i.e., halo sizes) for lensing (mostly 
early-type) galaxies. For this use of the outer power-law radial index, equation (1) will, in 
general, have a diverging total mass without a well-defined extent of the galaxy, although it 
would not be an issue in applications of the model to gravitational strong lensing, in which 
the extent of the dark matter halo is assumed to be much larger than the critical radius of 
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the lens. In order for the ‘cusp’ model [eq. (1)] with u 0 < 3 to have a well-defined extent 
and converging total mass, we could perhaps introduce an envelope declining steeper than 
u 0 = 3. In practice, it is likely that such a model is difficult to deal with mathematically. 
However, we could mimic an envelope using a subtraction term; in other words, we could 
use the following modified version of equation ( 1 ): 

p(R) = poR- Ui [(l + R 2 )-l v °- v *V 2 - (R 2 t + # 2 )-K-«'d/ 2 ]_ (39) 

For large radius, the above model [eq. (39)] can be approximated by p{R) oc R~^ 0+2 \ which 
permits a hnite total mass for any v a > 1 . 

Finally, it is worth noting that in this paper we limited ourselves to the families of two 
power-law mass distributions given by equations (1) and (39) mainly because of straightfor¬ 
ward mathematical tractabilities of the models. The following model [see Hernquist (1990) 
and Zhao (1996)], which is the generalized version of equation (39), would eventually be of 
use in future gravitational lens studies: 

p(R) = p 0 R~ Vi [(l + R?y^-^)h _ (^ + (40) 

where parameter 7 represents the “sharpness” of the break. However, until practical solutions 
to lensing by the model given by equation (40) are obtained and observational data sensitive 
enough to constrain parameter 7 necessitate use of the model, equation (39) corresponding 
to 7 = 2 in equation (40) will be useful to study galactic radial profiles by applying the 
model to gravitational lenses with the help of the calculational method presented in this 
paper. 

The author thanks Ian Browne, Chuck Keeton, Shude Mao, and Peter Wilkinson for 
encouragements, interests and useful comments on the manuscript. Chuck Keeton is also 
thanked for testing the code based on the results of this paper, which helped to find a 
typographic error in the code. The author thanks David Turnshek for the past and on-going 
encouragements and support. He also would like to thank Wyn Evans and Neal Jackson 
for interests and comments on this work. Financial support for this work comes from the 
universities of Manchester and Pittsburgh. 


A. Calculation of /-functions of ‘Nuker’ Model 

The /-functions of any surface mass density which is a function of ( [eq. (3)] are defined 
as follows: 

1 /* 27 t nr 

/ ( °)(r) = — / dcj)' / dr'r'niO, (Al) 

rn Jo Jo 
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where £ = £(r', 4>') = (A/r b )[l + ecos2(0' — 0 O )+ 2 and m = 1, 2, 3,.. 7 

Now for the ‘Nuker’ model [eq. (5)] for non-integer values of A we have 
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where g(0') = P + Q cos 2(0' — 0 O ) = r+[l + e cos 2(0' — 0 O )] (e, Q > 0) as defined in CKT. 
The integral in the last line of equation (A 6 ) can be evaluated, using equation (B3) of CKT, 
as follows 
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Now using equation (A4) of CKT and noting the definition of 4m( r ) above [eq. (A3)], 
equation (17) for i — 1 is obtained. 
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7 Notice that d[P°\r)\/dr = [— P°\r) + P° \r)\/r, and similarly for I^{r) and 4m( r )- 
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The hypergeometic function in equation (A8) can be re-written, via a linear transformation 
of hypergeometric functions (e.g. Gradshteyn & Ryzhik 1994), as follows 
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Substituting equation (A9) into equation (A8) and integrating over <f>', we find that the 
second term in equation (A9) vanishes. The rest of the calculation is very similar to that 
for -£> 2 m(r) and we obtain equation (17) for i = 2. The function /^ ) '( ker ' 3) (r; A, /i) [eq. (17) for 
i — 3] can also be obtained in a very similar way from the definition [eq. (A5)]. 


B. Revisiting ‘SPLE’ Model: Some New Expressions for /-Functions 

As described in §3, /-functions of the ‘cusp’ model [eqs. (2a,b)] are calculated using 
those of the ‘splc’ model [eq. (4)] as well as those of the ‘Nnker’ model [eq. (5)] since the 
‘cusp’ model can be expanded in terms of the ‘sple’ and ‘Nnker’ models. Here we consider 
some new expressions for the /-functions of the ‘sple’ model, which are useful for calculating 
/-functions of the ‘cusp’ model. 

The /-functions of the ‘sple’ model obtained in CKT converge very quickly either for 
large or small radii compared with the core radius, which is denoted by ry in equation (4). 
However, neither expression for the A 11 -function given in CKT converges sufficiently fast 
to be useful for calculating the /^-function of the ‘cusp’ model for r ~ ry, although the 
AT-function given in CKT is still useful. 

We now obtain new expressions for the A 1 )- and A 2 ^-functions of the ‘sple’ model which 
converge very quickly for r ~ r^. This is done via the following expansion of the ‘sple’ model, 

xc sple (r, 0) = (l + 
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After straightforward algebra using equation (Bl), we obtain 
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where f and Icos(/q m) are given in equation (30a). 


C. Prescriptions for Evaluating Non-trivial /-Functions 


All the /-functions obtained in CKT and in this paper are well-defined converging 
series. Numerical implementation of them is relatively straightforward. However, a few of 
them need to be dealt with carefully due to their non-trivially converging properties for 
part of the parameter space. Below we discuss these non-trivial /-functions and give some 
mathematical prescriptions which are useful for numerically evaluating them. 

Let us first consider the /^-function of the ‘splc’ model obtained in CKT, which is 
given by 


7 -sple( 2 ) 

1 2m 


('C k) 


ft(r)hV ^ rh(rr _ 

OO 

x y^[g 2 {r)] k F[m - k - n,-k - 1; £i(r)]. 

fc =0 


(Cl) 


While this function converges quickly for r > r b since £ 2 ( 7 ") — 1 ► 0 for r > r b , it converges 
slowly for r < r b especially as r/r b —> 0. The latter case (i.e. r < r b ) is irrelevant in most 
applications of the ‘sple’ model since observed images are normally found well outside any 
lens core radii. However, as can be seen in equations (14) and (23), evaluation of the / (2 h 
function of the ‘cusp’ model for r < ri, requires evaluation of the /^-function of the ‘splc’ 
model [eq. (Cl)] for r < r b multiple times. It is thus important to evaluate equation (Cl) 
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fast and accurately for r < r&. Fast and accurate evaluation of equation (Cl) can be done 
by expanding the hypergeomtric function and noting that infinite sum for each term in the 
hypergeometric function can be broken down to terms including 

OO 

S n (x) =Y J k n x k (n = 0,1,2,...), (C2) 

k =1 

which can be analytically evaluated using the following recursion relation 

— 1 n_1 / \ 

s„(z) = — £(- 1 )"-* ( I ) ( C3 > 

X k =0 k ' 

and S 0 (x) — x/(l — x). 

As noted in §3, the /(^-function given by equation (30a) of the ‘cusp’ model requires a 
special treatment. Equation (30a), which was obtained by substituting equation (B2) into 
equation (26), has a slow convergence over / although its convergence over k is fast. This 
slow convergence problem cannot be overcome using equation (27). In fact, the 1 ^-function 
given by equation (30a), which is needed for r ~ ry, is the only true slowly converging series 
out of all the /-functions given in CKT and in this paper. Let the Z-tli term in the infinite 
sum over l in equation (30a) be denoted by ui, i.e., 

u t = D t (ui, Vo )(1 - t )i r (^ + k + m + F + + k + m + n + 1,1; k + 2m + 2; t). (C4) 

T(l + /i + 1) 

As / —> oo, we have ui/ui+i —> 1 and 
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Table 1 


Summary of Mass Models 


name 

meaning 

projected density 

3-dimensional density 

SPLE 

softened power-law ellipsoid 

Equation (4) 

P(R) = P o(i + R 2 )-" /21 

Nuker 

“Nuker-law” model 

Equation (5) 

Not Applicable 

Cusp 

cusped two power-law model 

Equations (2a,b) 

Equation (1) 


1 See equation (1) for the definition of R. 







